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ABSTRACT 

The subject of this paper is the secular behaviour of a pair of planets evolving under 
dissipative forces. In particular, we investigate the case when dissipative forces affect 
the planetary semi-major axes and the planets move inward/outward the central star, 
in a process known as planet migration. To perform this investigation, we introduce 
fundamental concepts of conservative and dissipative dynamics of the three-body prob- 
lem. Based on these concepts, we develop a qualitative model of the secular evolution 
of the migrating planetary pair. Our approach is based on analysis of the energy and 
the orbital angular momentum exchange between the two-planet system and an exter- 
nal medium; thus no specific kind of dissipative forces is invoked. We show that, under 
assumption that dissipation is weak and slow, the evolutionary routes of the migrating 
planets are traced by the Mode I and Mode II stationary solutions of the conservative 
secular problem. The ultimate convergence and the evolution of the system along one 
of these secular modes of motion is determined uniquely by the condition that the 
dissipation rate is sufficiently smaller than the proper secular frequency of the system. 
We show that it is possible to reassemble the starting configurations and migration 
history of the systems on the basis of their final states and consequently to constrain 
the parameters of the physical processes involved. 
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1 INTRODUCTION 

The secular regime of motion of multi-planetary systems is 
universal; in contrast with the 'accidental' resonant motion, 
characteristic only for specific configurations of the planets, 
secular motion is present everywhere in phase space, even 
inside the resonant region. The secular behaviour of a pair 
of planets evolving under dissipative forces is the principal 
subject of this paper, particularly, the case when the dissi- 
pative forces affect the planetary semi-major axes and the 
planets move inward/outward the central star, the process 
known as planet migration. It is presently accepted that mi- 
gration plays an important role in the dynamical history of 
planetary systems. Several mechanisms have been proposed 
to be responsible for planet migration; some of these are i) 
planet interactions with a protoplanetary disk of gas/dust, 
ii) gravitational scattering and clearing of remnant planetes- 
imal debris by the planets, iii) direct collisions between the 
planets, and iv) tidal interactions of the planets with a cen- 
tral star. There exists a vast literature on this issue; the 
reader is referred to Armitage (2010) for reviews of planet 
migration and references therein. 

The study of planet migration is frequently connected 
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to theories of planet formation. The modeling of interac- 
tions of growing planets with protoplanetary discs is a high- 
complexity task independently whether it is done analyti- 
cally or in the form of numerical simulations. Despite the 
significative progress stimulated by the discovery of extra- 
solar planets, the studies still suffer from various limitations 
(see Lubow and Ida (2010) for the recent review). The main 
obstacle is probably a lack of knowledge about the structure 
of planet-forming gaseous discs and their detailed physical 
properties. The huge list of unknown parameters requests 
numerous numerical simulations, characterized by high com- 
putational costs, and introduces uncertainties in the results 
obtained. Together with poorly determined physical proper- 
ties of the process, unknown starting configurations of the 
planets complicate investigations in the cases of migration 
of the tidally affected close-in planets and scattering the 
planetesimals in the late stage of the planetary formation. 

In this paper, we purpose a simple method for a quali- 
tative study of the planet migration, which does not require 
a detailed knowledge of dissipative mechanism, its physical 
properties and starting configuration of the migrating sys- 
tem. This is because we do not consider any specific kind of 
dissipative forces (they can be originated by tidal torques in 
the disc, tidal interactions with the star, ejection of planetes- 
imals by planets, direct collisions etc). Knowing that, under 
dissipation, the energy and the orbital angular momentum 
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of the planet system are no longer conserved, we model the 
migration through the variation of the energy and parame- 
ters of the secular dynamics: the angular momentum and the 
semi-major axes ratio. Similar approaches have already been 
used successfully in dynamical studies (see Malhotra 1995, 
Murray et al. 2002, Nelson and Papaloizou 2002, among oth- 
ers). 

Exploring basic concepts of the conservative and dissi- 
pative dynamics, we develop a model of secular dynamics 
of a migrating pair of planet in Section [2] We first ana- 
lyze possible effects of dissipation on the global dynamical 
quantities of the secular system, such as energy and orbital 
angular momentum, which are invariable during the con- 
servative motion. From these effects, the most known is the 
increase/decrease of the orbital energy of the system defined 
by the Keplerian motions of the planets. The result is the 
variation of the semi-major axes and the consequent expan- 
sion/contraction of the planetary orbits. 

The effects of dissipative forces on the secular compo- 
nent of the total energy of the system have been studied 
much less. The secular energy is defined by mutual secular 
interactions of the planets and its dissipation affects mainly 
the planet eccentricities provoking their damping oscilla- 
tions. Even when only one planet is directly subjected to 
friction forces, the orbital angular momentum exchange be- 
tween the planets affects the eccentricity of the other planet. 
We show that, if the dissipation is sufficiently weak and 
slow (adiabatical approach), the planets ultimately evolve 
into a nearly steady state defined by a stationary solution 
of the conservative problem, either Mode I or Mode II of 
motion (Michtchenko and Ferraz-Mello 2001), depending on 
the physical parameters of the system. In what follows, the 
already damped system continues to evolve under acting dis- 
sipative forces, tracing the family formed by the stationary 
solutions. We show that the final configuration of the system 
is defined by the balance between the orbital and secular en- 
ergy variations during migration, together with the angular 
momentum exchange between the system and an external 
medium. 

The developing of the model of the migrating pair of the 
planets follows essentially two steps: The first one is the cal- 
culation of the families composed of the Mode I and Mode 
II stationary solutions, for all possible values of planetary 
semi-major axes and angular momentum of the system. Pa- 
rameterized by the mass ratio, these families will define the 
evolutionary routes of the migrating system in phase space. 
To obtain the families of stationary solutions, with no re- 
strictions on the values of planet eccentricities, we use the 
precise semi-analytical approach developed in (Michtchenko 
and Malhotra 2004). 

The second step consists of the investigation of the sta- 
bility of Mode I and Mode II of motion, to determine which 
of these modes will play a role of the center, toward which 
the system is attracted during migration. We show that the 
stability of a center depends on the combined effect of three 
conditions: The first one is related to the law of the orbital 
angular momentum exchange between the system and the 
external medium. For instance, it can imply that the orbital 
decay of the planets results in the decreasing of the eccen- 
tricities of their orbits and vice versa. The second condition 
is related to the loss/gain of the secular energy by the sys- 
tem. It should be noting that the gain of the secular energy 



results in divergent planetary orbits, while the loss of the 
secular energy in convergent orbits. The third condition re- 
lates two parameters of the secular problem, the mass ratio 
mj/mi and the semi-major axes ratio ai/a2 (indices 1 and 2 
correspond to the inner and outer planets, respectively). We 
show that Mode I of motion plays a role of an attractive cen- 
ter, when planetary orbits diverge and Wai/a? < rrn/mi, 
or, when planetary orbits converge and \fa\fa% > ms/mi. 
In the opposite cases, the Mode II of motion is stable. We 
show that an immediate consequence of this result is that 
the evolutionary history of a dissipative system can not be 
accessed through the back in time simulations of its evolu- 
tion. 

In order to illustrate the performance of our model, we 
do a series of purely numerical experiments, whose results 
are presented in Section [3] We apply the model to several 
fictitious systems, varying the planetary masses and semi- 
major axes and simulating both types of orbits, divergent 
and convergent ones. We initially assume that only the in- 
ner planet is affected by external dissipative forces and that 
the orbital angular momentum of the system is conserved 
during migration. It should be noted that these assumptions 
describe well the secular interaction of a tidally affected syn- 
chronously rotating planet with its outer companion. 

We test then the dependence of the migration evolution 
on the starting conditions of the secular system. We verify 
that the final configuration of the system is independent 
on where in the phase space the system started its secular 
evolution, provided that the starting values of eccentricities 
satisfy the conservation of the orbital angular momentum. 
Thus, despite a qualitative nature of our model, it is able 
to indicate precisely the domains in phase space where the 
system could start its dissipative evolution. We also inves- 
tigate the dependence of the simulated planetary paths on 
the rate of the migration process and the individual masses 
of the planets and determine the applicability limits of the 
adiabatical approach. 

In the next section, we extend the model to the case 
when the system exchanges orbital angular momentum with 
the external medium. We introduce the AM-leakage, de- 
fined as a portion (a) of the angular momentum variation 
produced by migration (i.e. the expansion/cotraction of the 
planetary orbits), which is extracted (or added to) from the 
system. We construct the migration routes of the system 
evolving with different values of the factor a. This approach 
provide us with a general idea of how the system could evolve 
under a variety of migration conditions and physical models. 

In Section [5] the model is adapted to the case when 
only the motion of the outer planet is affected by dissipa- 
tive forces. Analyzing the migration tracks obtained in this 
case, together with the previously obtained ones, we show 
that the evolution and the final state of the secular system, 
evolving under weak dissipative forces, with typical values of 
q smaller than 1, depend on the balance between the orbital 
and secular energy variations, as follows: 

- the loss of the orbital and secular energy defines the 
total circularization of two convergent orbits, allowing a cap- 
ture in low-order mean-motion resonances; 

- the loss of the orbital energy and the gain of the secular 
energy define the totally circularized orbits and the increas- 
ing mutual planetary distance; 



Secular model of migrating planet pairs 3 



- the gain of the orbital and secular energy defines the 
increasing mutual planetary distance, together with the in- 
crease of the eccentricities of the orbits; 

- the gain of the orbital energy and the loss of the secu- 
lar energy define the continuous increase of the eccentricities 
of two convergent orbits, when the secular system is broken 
due some external perturbations, such as mean-motion res- 
onance interactions or close encounters between the planets. 

In Discussions, we show several examples found among 
the planets of the Solar System and the known extra-solar 
systems, whose actual configurations could contain records 
of the dissipative evolution in the past, according to one 
of the above scenarios. As a final consideration, we discuss 
advantages and disadvantages of the model developed in this 
paper. 

Fundamentals needed to understand the developed 
model are given in Appendices [X] and [B] In Appendix [X] 
we discuss the main features of the conservative secular mo- 
tion, developing the first-order planar Hamiltonian model of 
the three-body problem. We show that: i) it is linear when 
the motion is projected on a sphere (Pauwels 1983); ii) the 
planet mass and semi-major axes ratios, together with the 
angular momentum of the system, play the role of free pa- 
rameters of the secular problem; iii) the eccentricities of the 
planetary orbits oscillate in anti-phase around one of two 
stable stationary solutions, known as Mode I and Mode II of 
motion. Recall that Mode I characterizes the aligned config- 
uration of the planetary orbits, when the averaged mutual 
distance between the planets is minimal and consequently 
the energy of the system is also minimal. On contrary, Mode 
II characterizes the anti-aligned configuration and the max- 
imal energy of the system, when the averaged mutual dis- 
tance between the planets is maximal. 

In Appendix [B] we study the behaviour of a linear sec- 
ular system moving under action of dissipative forces. We 
show that the behaviour of the system is similar to the be- 
haviour of the well-studied harmonic oscillator moving under 
friction (see Andronov et al. 1966, Chap. I). The anti-phase 
oscillations of the planet eccentricities around a fixed center 
are damped and the system ultimately evolves into a steady 
state defined by a stable focus. We show that the stable fo- 
cus of the secular system is close to the center defined either 
Mode I or Mode II stationary solutions of the conservative 
problem, depending on the physical parameters of the sys- 
tem and dissipation. We note that this result is independent 
on the specific form of the dissipative force, but requires only 
sufficiently small and slow dissipation of the energy. 



2 THE MODEL 

2.1 The conservative secular evolution 

First we consider a conservative three-body system consist- 
ing of the central star (M) and two planets with masses mi 
and ni2 (hereafter, the indexes i = 1, 2 correspond to the 
inner and outer planets, respectively). In the heliocentric 
reference frame, the canonical set of variables is constituted 
by the relative position vectors fi and conjugate momenta 
Pi = m,i where the position vectors pi are relative to 
the center of gravity of the three-body system. If only mu- 
tual gravitational interactions are taking into account, the 
Hamiltonian of the problem is written in the form 
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where k 2 is Gaussian constant, fii — k 2 (M + rrii), m'i — 
rriiM/(M + rm) and A = |ri — 7*2 1 . The first term is the 
sum of Keplerian motions of the planets and the second and 
third terms produce direct and indirect perturbations among 
the planets, respectively. 

We suppose that the planets move in the same plane; 
then the Hamiltonian JTJ can be written, in terms of canon- 
ical elliptic variables and the disturbing function R, as 
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where the mass-weighted Poincare variables are introduced 

as 

Xi mean longitude, Li = 

— vji longitude of perihelion, Ii = 

at and ei stand, respectively, for the semi-major axes and the 
eccentricities of the planetary orbits. In order to investigate 
the secular behaviour of the system, we apply the averaging 
procedure to the Hamiltonian ([2]). The averaging done with 
respect to the mean longitudes of the planets, removes from 
the problem the short periodic oscillations. The averaged 
Hamiltonian is defined by 
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where the averaging procedure need be done over only the 
direct part of the disturbing function, because the indirect 
part does not contain secular terms (Brouwer and Clemence, 
1961). 

After the elimination of the short periodic terms, the 
averaged Hamiltonian does not depend on Xi; consequently, 
the semimajor axes of the planet orbits, ai and a%, are 
constant and serve simply as parameters in the Hamilto- 
nian ©. In addition, due to the D'Alembert's rule, the 
^-dependence in the averaged disturbing function is only 
through At37 = tU2 — zo\ ; this implies that I\ + I2 is also a 
constant of motion. This quantity is so-called angular mo- 
mentum deficit of the system and, up to second order in 
masses, can be written as 



AMD = y^rrtj ■ ; 



(5) 



In contrast with the orbital angular momentum of the sys- 
tem given by 



AM 



r t a 2 a/T~ 



(6) 



AMD is invariable solely in the secular problem, due to the 
fact that, in this problem, a\ and 0,2 are parameters. 

It is worth emphasizing that, since the secular inter- 
actions do not depend on the positions of the planets on 
the orbits, the secular behaviour of the system can be de- 
scribed by gravitational interactions between two coplanar 
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rings, each one with the mass of the corresponding planet 
(Wu and Goldreich 2002). The massive rings do not change 
their sizes under secular perturbations (at = const), but only 
their shapes, defined by a, and the mutual orientation, de- 
fined by the secular angle Azu. The detailed description of 
the conservative secular dynamics of the two-planet system 
is given in Appendix 1X1 

The energy of the system given by the Hamiltonian Q 
is conserved during the secular evolution. For the purpose of 
this work, it can be re- written in analogy with the expression 
CQ) as: 



— _ _fnm 1 



2ai 
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2a 2 
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where the indirect part of the disturbing function was omit- 
ted. In the above equation, the two Keplerian terms are 
dominant; they are functions solely of the parameters ca 
and, consequently, are constant in time. Thus, due to the 
conservation of the energy of the system, the latter term 
must be also conserved during the evolution. The constant 
D is introduced as 



D 



-k mim 2 /ri a 



(8) 



which is an increasing function of the increasing semi-major 
axes ratio a\ja% (Callegari et al. 2004). 



2.2 The dissipative secular evolution 

Now we suppose that two planets, interacting secularly (the 
planets are far enough from any mean-motion resonance), 
undergo external perturbations which affect the energy and 
angular momentum of the system. One consequence of such 
perturbations is that the planet semi-major axes are altered 
and the planets move inward/outward the central star, the 
process known as planet migration. The variation of the en- 
ergy ([7]) over a small time increment is given by 



aT7 fiim'iAtti ^2JTi2^ a 2 k 2 m\m 2 AD 
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where Aai, Aa 2 and AD are variations of the semi-major 
axes and the separation between the orbits, respectively. The 
sum of the two first terms defines the change of the orbital 
(Keplerian) energy of the system, while the latter term de- 
fines the change of the secular energy of the system. In what 
follows, we will show that the balance between both defines 
the migration path of the system in phase space. It is worth 
noting that the variation of the secular energy can occur 
even when the semi-major axes are not affected by external 
forces, but when eccentricities are changed (see Appendix 

ED. 

For the sake of simplicity, we suppose initially that only 
the semi-major axis of the inner planet a± is affected. The 
orbital decay of the close-in inner planet due to tidal inter- 
actions with the central star can be approximated by this 



where % scc is the secular part of the Hamiltonian given in 
Equation Q, in general form, and in Equation (|AlJl . in the 
linear approximation. Despite the fact that D is a compli- 
cated function of the planet semi-major axes and eccentrici- 
ties, it can be used as a convenient measure of the separation 
between two orbits. For instance, in the case of two circular 
orbits, D = a2/by 2 , where &J/ 2 is the Laplace coefficient 



assumption. (The case of the change of a 2 will be discussed 
in Section [5]) 

According to Equation @, the loss of orbital energy 
of the inner planet (Aai < 0) does not affect the orbital 
energy of its companion; thus, the semi-major axis of the 
outer planet remains constant during migration. In the other 
words, there is no exchange of the orbital energy between 
the planets in the secular approximation (Wu and Goldreich 
2002). The coupling between two orbits is defined by the 
secular term in Equation JSJ; it is responsible for transfer 
of the tidally induced changes in the orbital elements of the 
inner planet to the outer planet (for details see Appendix 
B2|) . It is worth noting that, during divergent migration, 
when Aai < (or Aa 2 > 0), the system gains the secular 
energy, since the separation between two orbits is increased 
(AD > 0). Inversely, during convergent migration, when 
Aai > (or Aa2 < 0), AD < and the system loses the 
secular energy. 

We further assume that the orbital angular momentum 
of the system is conserved during migration. In the approxi- 
mation that only the inner planet is affected by tidal interac- 
tions with the star, this assumption is valid when the planet 
rotation has been synchronized and effects of the stellar tides 
are not considered (the case of tidal interactions of a close-in 
planet with a very slow rotating central star) (Casenave et 
al. 1980, Correia and Laskar 2010, Rodriguez et al. 2011a). 
In this case, the decrease Aai produces the damping of the 
inner planet eccentricity as follows: 



(1 



2aiei 



^Aai, 



(10) 



where the index "ex" indicates that this eccentricity varia- 
tion is a response to external forces and it must be added 
to the variation originated by the secular interaction with 
the other planet. Note that according to the above equation, 
under the assumption of invariable AM, the orbital decay is 
halted when the planet orbit is circularized. It should be em- 
phasized that the hypothesis of the invariable orbital angular 
momentum is introduced here to facilitate our understand- 
ing of the dissipative behaviour of the system; in Section 
[4] the model will be extended to more general case, when 
the angular momentum of the system is exchanged with the 
external medium. 

In contrast with the orbital energy, there exists the ex- 
change of the orbital angular momentum between the plan- 
ets. From Equation ©, we can express, up to second order 
in masses, the coupling variations of both eccentricities as 



Ae 2 




iAei. 



(11) 



The above equation shows that, in addition to the secular 
anti-phase variations of both eccentricities, the orbital de- 
cay of the inner planet (Aai) and the eccentricity damping 
Aef x (|10p provoke the smooth decrease of e 2 (see for details 
Rodriguez et al. 2011b). 

The secular behaviour of the eccentricities under act- 
ing friction forces is analyzed in detail in Appendix [B] We 
show that during the orbital decay of the inner planet, both 
eccentricities present damped oscillations; when the ampli- 
tude of the oscillation of e\ tends to zero, the amplitude of 
e 2 also tends to zero. The only condition needed is that the 
dissipation rate is sufficiently slow, compared to the proper 
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Figure 2. Domains of the stable centers defined by the condition 
( l 1 2 I t on the parametric (n2/ni,rrt2 /mi)— plane, where n is the 
mean motion of the planet. This case corresponds to the condition 
(113ft ; in the opposite case, the stability of the domains is inverted. 
Locations of the main mean-motion resonances are indicated. 



Figure 1. The families of stationary solutions (Mode I and 
Mode II), parameterized by the mass ratio m2/m\ (values be- 
side each curve). Top: the (ei,e2) representative plane. Bottom: 
the (ei,rai/ri2)-plane. On both panels, Aro is fixed at along pos- 
itive values on the ei-axis) or 180° along negative values on the 
ei-axis). The branches of the curves in black color corresponds 
to solutions with if < I,f , and in red color to solutions with 
if > 1 2 , where if and iS are the partial angular momentum 
deficits of the planets. Arrows show the direction of divergent 
orbits. 



period of the secular oscillation. The final configuration of 
the damped system is either aligned or anti-aligned orbits, 
depending on the sign of Ae"" and the physical parame- 
ters of the system, the mass and semi-major axes ratios. 
Seen in phase space (e.g., Figure iBl]) . the damped system 
moves along the spiral trajectory attracted to the stable cen- 
ter, which can be either Mode I (aligned orbits) or Mode II 
(anti-aligned orbits) stationary solution of the conservative 
problem, characterized by the maximal and minimal energy, 
respectively. This state of the system is often referred to as 
'quasi-stationary' in the literature (Mardling 2007, Batygin 
et al. 2009). 

During the evolution, the attracting center itself dislo- 
cates slowly in the phase space, due to changes of the secular 
parameter axja-z and the consequent change of the orbits 
separation D (see an example in Figure IB3[) . Even when 
'quasi-stationary' state is attained, the system continues to 
evolve under action of dissipative forces, following closely 
the Mode I or Mode II family of stationary solutions of the 
conservative problem (some of these solutions are shown in 
Figure [T]). The smooth migration of the system is halted 
only when either both orbits are circularized, in the case 
of Aoi < 0, or the system enters into the domain of phase 
space where secular model is not longer valid, in the case of 
Aai > 0. In the former case, Aei x = Aai — 0, when ej = 0, 
according to Equation (|10[) . In the latter case, the planets 
approach each other sufficiently to introduce strong short- 
terms or resonance perturbations which are not considered 
by the secular theory. 



2.3 Evolutionary routes 

As described above, the families of the conservative station- 
ary solutions draw migration routes of the secular systems 
evolving under friction forces. Three of those families are 
shown in Figure [T] The families were calculated using the 
algorithm described in Michtchenko and Malhotra (2004), 
for continuous values of the ratio ai/a<2 and for three dif- 
ferent values of the mass ratio m^/mi, namely 2.0, 0.5 and 
0.3. Along each family AM is fixed at the value calculated 
for the system given by a\ = 0.02 AU, = 0.1 AU and 
ei = e 2 = 0.01. 

We plot the families on two different planes: the (ei,e2)- 
plane of the planet eccentricities (top panel) and the 
(ei,ni/ri2)-plane (bottom panel), where the ratio of the 
mean motions of the planets ni/nz is computed from the 
ratio of their semi-major axes through the third law of Ke- 
pler. On both panels, Aw is fixed: at 0, along positive values 
on the ei-axis, or at 180°, along negative values on the ei- 
axis. 

Each family consists of two distinct branches connected 
at the origin. The branch with positive ei-values corresponds 
to Mode I of motion, while the branch with negative ei- 
values corresponds to Mode II. Which from these modes will 
play a role of the stable center and, consequently, trace the 
trajectory of the migrating system is discussed in detail in 
Section [B3I As shown, the stability of the center is defined 
by the relation between if and 1% , which are the partial 
angular momentum deficits of the planets, evaluated at the 
center. For Aai < (and Ae™ < 0), the system evolves 
into the center with if < 1% . This is also true for Aa2 > 
(and Ae| x > 0). In both cases (orbital decay of the inner 
planet and orbital expansion of the outer planet) the planets 
move away from each other. The orbits are diverging, their 
separation D is increasing, and the system gains the secular 
energy, according to Equation ©. 

The situation is inverse for convergent orbits, when 
Aai > or Aa2 < 0. In this case, the separation D is 
decreasing and the system loses the secular energy. As a 
consequence, it evolves into the center with if > I 2 ■ The 
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transition condition I{ = 12 
order in masses and eccentricities, as 



I§ can be re-written, up to second 




(12) 



In Figure [T] the centers with if < if are plotted 
by black dots and, those with if > if, by red dots. Ac- 
cording to described above, for divergent orbits, the Mode 

I of motion is stable for ma/mi = 2.0, while the Mode 

II is stable for rri2/mi = 0.3. The transition value of the 
mass ratio given by the condition (|12[) is equal to 0.447, for 
0,1/0,2 = 0.2. The family of stationary solutions parameter- 
ized by rm/mi = 0.5 in Figure [TJ is close to the transition 
case. We can observe that, in this case, the centers change 
the stability during the evolution of the system along the 
branch. Thus, in order to access the stability of the cen- 
ter, when 7712/mi is close to the critical value, the condition 
(|A9|) . which takes into account planetary eccentricities, must 
be applied. 

It is worth emphasizing that the determination of sta- 
bility of the center was done assuming that Aai and Aei 
have the same sign (see Equation (|10p ). that is a consequence 
of the conservation of the orbital angular momentum of the 
system during migration. In the more general case, we in- 
troduce the condition 



sgn(Aai, 2 ) = sgn(Aei 



(13) 



and resume the information on the stability of the centers in 
Figure [2] where we plot the curve given by Equation (|12[) . 
on the parametric plane {nz/ni, mz/mi). If the condition 
(| 13p is satisfied, the divergent orbits ultimately converge to 
the Mode I stationary solutions in the domain above the 
curve. In the same domain, the convergent orbits tend to 
the Mode II solutions. In the domain below the character- 
istic curve, the stability of the secular modes is opposite. 
If the condition (| 13f) is not satisfied, the stability shown in 
Figure [5] is just inverted. It should be noted that the con- 
dition (|13l) is typical for the known dissipative processes in 
the planetary systems. The dynamical interpretation of this 
condition, from the point of view of the orbital angular mo- 
mentum exchange, will be discussed in Section [4] 

According to our model, the system moving under slow 
dissipation evolves following the stationary solutions of the 
conservative secular system, which draw the evolutionary 
route that the planets acquired from their initial configu- 
rations towards the present locations. Three such families 
parameterized by different mass ratios were shown in Fig- 
ure [T] Analyzing these routes and migration conditions, we 
can i) do robust predictions about where in the phase space 
the system could arrive from, and ii) conjecture on the pos- 
sible final configurations of the system. For instance, in the 
case of inner planet migration, we conclude that, when the 
system loses its orbital energy and gains the secular one, its 
evolution is characterized by the divergent migration and 
decreasing eccentricities, with final configurations on the to- 
tally circularized orbits. The value of a\ has a lower limit 
which can be obtained from Equation (J5]), for given AM and 
02 and fixing ei = t2 — 0; it is given in Equation (|B3[) . 

During the convergent migration, when system gains 
the orbital energy and loses the secular one, the evolu- 
tion is characterized by increasing eccentricities; the final 
state of the system is uncertain, since it approaches the 
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Figure 3. Dependence of the migrating evolution on initial con- 
ditions on the (ex 1 e2)-plane. Three different starting positions 
of the system are shown by cyan symbols. The final (current) 
configuration of the system is shown by a black symbol. Families 
of stationary solutions are plotted by thick red curves. Mode I 
corresponds to positive ex— values, while Mode II to negative ex- 
values. Arrows show the direction of the evolution of the system. 
The mass ratio is fixed at mz/mi = 2.0. 



strong mean-motion resonances and close encounters do- 
main, where short-period perturbations could disrupt the 
system. It should be emphasized that the described scenar- 
ios are possible when the migration condition (|13|l is satis- 
fied. In the opposite case, the behaviour of the eccentricities 
is inverted. Finally, the similar analysis with respect of the 
migrating outer planet will be done in Section [S] 



3 NUMERICAL SIMULATIONS 

To illustrate the behaviour of the two-planet secular system 
evolving under dissipative forces, we perform a series of nu- 
merical simulations. For this task, we use a purely numerical 
approach which consists of integrations of exact equations of 
motion averaged through the applied on-line low-pass filter 
(Michtchenko et al. 2002). We do not invoke any specific kind 
of dissipative forces; the orbital migration is modeled in the 
exact equations of motion of the planets through constant 
perturbations in their semi-major axes and eccentricities, in 
analogy with Malhotra (1994) and Nelson and Papaloizou 
(2002). 

The system under study consists of the central star 
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with mass M = 1.0 Ms un and two gravitationally interacting 
planets. Different values of the planetary masses were used in 
this work; they will be specified for each experiment reported 
in this section. We place the system in its current configu- 
ration described by the semi-major axes ai = 0.02 AU and 
a,2 — 0.1 AU and small eccentricities e\ — e2 = 0.01. In this 
section we assume that the dissipative forces are acting only 
on the inner planet. We also assume that the orbital angu- 
lar momentum of the system is conserved (the more general 
case of the AM transfer will be studied in Section [4| and 
calculate its value for the current configuration of the sys- 
tem. Thus, during the simulation, the inner semi-major axis 
is forced to decay by Aai, after each time increment (typi- 
cally smaller than one tenth of the orbital period), and the 
value of ei is corrected by the amount given in (|10p . For the 
sake of simplicity, in this work we adopt a linear regime of 
the orbital decay, in the form < di >= const, where const 
may be positive or negative. 

It should be emphasizing that the model allows us to 
study separately the secular regime of motion of the sys- 
tem. In fact, the averaging procedure shown in Equation 

removes all effects of short-period terms and of mean- 
motion commensurabilities. This is not a case of the nu- 
merical integrations of the exact equations of motion, whose 
output contains information on all possible regimes of mo- 
tion, which are overlapping intrinsically. This fact should be 
kept in mind when the results obtained by both approaches 
are compared. 

3.1 Dependence on initial conditions 

We first test the dependence of the migration evolution of 
the secular system on its initial configuration. We fix the set 
of the planetary masses at mi = 1.0 Mj and m.2 = 2.0 Mj, 
and choose arbitrarily the starting semi-major axes as a± = 
0.0311 AU and a 2 = 0.1 AU. These values correspond to 
n\/ri2 = 5.76 and place the system far from any strong 
mean-motion resonance. 

The description of the starting configuration of the sys- 
tem must be completed by the initial values of ei, e2 and 
Acer. Their choice is illustrated in Figure [3l where we show 
three representative planes (ei,e2), each one presenting the 
different starting position of the system. Since the orbital 
angular momentum is conserved in this experiment, all pos- 
sible initial values of eccentricities are located on the AMD- 
level obtained, for given planetary masses and initial a\ and 
(Z2, from Equation ((5J; this level is shown on all graphs in 
Figure together with the families of stationary solutions 
plotted by thick red curves. 

The first set of the eccentricities (cyan symbol on the 
top graph) was chosen in the neighborhood of the Mode 

1 solution, with the initial Aw = 0. The second set (cyan 
symbol on the middle graph) was chosen in the opposite con- 
figuration, in the vicinity of the Mode II solution, with the 
initial Aw = 180°. The third set of the initial eccentricities 
(cyan symbol on the bottom graph) was in the intermedi- 
ate configuration, with the initial value of Aw within the 
interval from to 360° . 

The chosen initial conditions were numerically propa- 
gated in time, according to the described above procedure, 
with the constant dissipation rate < di >= — 10 _9 AU/yr; 
the resulting trajectories are plotted in the (ei,e2)-planes 



by solid curves in Figure [3] In the box, on each panel, we 
show the time evolution during migration of the secular an- 
gle Aw. 

Figure [3] shows that, independently, where the system 
starts the evolution, the divergent orbits, with Aai < 0, ul- 
timately evolve into the Mode I of motion (Aw = 0), which 
guides the system towards its current position in the vicin- 
ity of the origin. This result, obtain ed for t he system whose 
parameters satisfy the condition -^/ a i/ a 2 < mi/m\, is in 
complete accordance with our model described in Section [5] 

The important conclusion from the performed experi- 
ment is that the final state of the migrating secular system 
is independent of its initial position in the phase space, de- 
fined by given values of the parameters: the mass and semi- 
major axes ratios and the orbital angular momentum. As a 
consequence, to assess the system configuration in the past, 
we have to know solely the starting value of ai, while the 
values of both eccentricities and of the secular angle may 
be chosen arbitrarily, with only restriction that the eccen- 
tricities belong to the AMD-level defined by the masses and 
the semi-major axes of the planets. Inversely, the obtained 
in such a way magnitudes of the eccentricities could impose 
constrains on the possible starting values of oi. 

3.2 Dependence on the mass and semi-major axes 
ratios 

In this section we test the stability of the domains shown 
in Figure [2] For this task, we vary the masses and semi- 
major axes of the planets and consider both types of mi- 
gration, divergent and convergent. The output of the four 
performed runs is shown in Figure [4] on the representative 
planes (ei,e2) (top row) and (ei,ni/ri2) (middle row). On 
each plane we plot the family of stationary solutions ob- 
tained for the planet mass ratio used in the simulation; the 
mass ratio values are shown on the top of each column. The 
time evolutions of the secular angle Aw are shown on the 
bottom panels, each one indicating the number of the cor- 
responding run. 

3.2.1 The case of \/oi7<i2 < ma /mi 

The first run (#1) simulates the migration of the system 
composed of two planets with the masses mi = 1.0 Mj 
and m,2 = 2.0 Mj. The starting inner semi-major axis was 
ai = 0.0455 AU, which corresponds to ni/ng = 3.24. Since 
the final configuration of the system is independent on the 
initial eccentricities, their values were chosen arbitrarily on 
the initial AMD-level, in the vicinity of the Mode I solution, 
with initial Aw — 0. The dissipation was simulated with the 
constant rate < di >= — 10~ 9 AU/yr, which induces the di- 
vergent orbits and the gain of the secular energy. 

The resulting path is shown by black dots on the graphs 
from the left column in Figure [4] Its starting position shown 
by a cyan symbol is indicated by #1. For given values of 
the masses and semi-major axes, the condition -^/ a i/ a 2 < 
1712/ mi is always satisfied during the evolution of the sys- 
tem. In this case, according to our model, the dissipation of 
the energy causes the damped eccentricity oscillations and 
the convergence of the system to the Mode I of motion. 
Consequently, this secular mode drives the system toward 
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Figure 4. Numerical simulations of the migrating two-planet system. Families of stationary solutions are plotted by thick red curves. 
Mode I corresponds to positive ei— values, while Mode II to negative ei— values. Starting positions of the system are shown by cyan 
symbols and indicated by the corresponding run number. Runs #1 and #4 are shown by dots, while #2 and #3 are shown by solid lines. 
Arrows beside each path show the direction of the evolution of the system. 



the origin, where both orbits are circularized. This scenario 
is fully reproduced by the simulation #1. The time evolu- 
tion of the secular angle Aw shown on the bottom panel 
in Figure 21 clearly illustrates the capture in the Mode I of 
motion. 

The observable deviations of the planetary path from 
the Mode I of motion is due to passages through sev- 
eral mean-motion resonances during migration. During these 
passages the planet eccentricities are strongly excited, but 
they are damped again as soon as the system leaves the 
mean-motion resonance. The correlation between the pas- 
sages through the mean-motion resonances and excita- 
tions of the eccentricities can be clearly observed on the 
(ei,ni/ri2)-plane in Figure ^middle. Far from the main 
mean-motion resonances, only the component of the smooth 
decrease can be observed in the evolution of the planetary 
eccentricities. 

It is worth noting that, in the theory of oscillations, the 
capture into the quasi-steady state is often referred to as 're- 
laxation', to indicate the return of the system to equilibrium 
after being perturbed by an external force. In this work, we 
do not use this term, because there is no qualitative differ- 
ence between finite-amplitude oscillations around the stable 



center and nearly zero-amplitude oscillations (and smooth 
evolution) in the very close vicinity of the center. 

Now we suppose that the position marked by #2 on the 
top and middle panels in Figure \$\left column, is the current 
position of the system. This configuration can be used as the 
starting point to simulate the migration of the same pair of 
planets, evolving now under the positive variation of the 
inner semi-major axis, with the rate < di >= 10 _9 AU/yr. 
When the orbit of the inner planet is expanded, the evolution 
of the system is convergent and the secular energy is lost. 
Note that the run #2 may be considered as a back in time 
simulation of the run #1. 

The trajectory obtained is shown by a continuous line 
on the graphs in Figure [4]Ze/i column. Starting on the nearly 
circular orbits with aligned periastra, the planet eccentric- 
ities are exited very quickly. The increasing amplitudes of 
oscillations lead to the circulation of the secular angle Aro, 
whose behaviour is shown on the bottom panel in Figure U 
After a stage of circulatory evolution, the periastra of the 
orbits enter in the anti-aligned regime of motion, when the 
system converges to the Mode II of motion. The system con- 
tinues to evolve along this secular mode in the direction of 
increasing eccentricities, crossing several mean-motion res- 
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Figure 5. Dissipativc evolution of the system along Mode I family 
of stationary solutions, for three different values of the dissipa- 
tion rate. The mass ratio 7712 /mi is fixed at 3 and all orbits start 
at the same initial configuration, shown by a black symbol. The 
excitation due to the passages through main mean-motion res- 
onances and damped oscillations around Mode I can be clearly 
observed. 



onances, until the convergent orbits reach the strong 4/1 
resonance. The chaotic effects of this resonance on the high- 
eccentricity planet motion cause the disruption of the sys- 
tem. 

We verify that the result obtained is in perfect accor- 
dance with our model, which pre dicts t hat the convergent 
orbits, satisfying the condition -J a\ Ja2 < rn^/mi, will ul- 
timately converge to and follow the Mode II of motion (see 
Figure [2| . It is worth noting that the evolution of the tra- 
jectory #2 is completely distinct from that of the run #1, 
despite that the run #2 is just a back in time simulation 
of the run #T. Thus, the important conclusion from the 
done experiment is that, for dissipative systems, forward and 
backward in time integrations are qualitatively different; in 
the other words, the evolutionary history of dissipative sys- 
tems can not to be studied through integrations going back 
in time. 



3.2.2 The case of ^Jai/a.2 > mzjmi 

According to our model illustrated in Figure [21 the stability 
of the centers shown in the previous section is inverted in 
this case. By definition, a\ < 02 always; thus the above 
condition is satisfied only for the systems with a less massive 
outer companion. For this reason, two simulations, #3 and 
#4, shown on the graphs from the right column in Figure SI 
set the planetary masses at m,\ — 1.0 Mj and 7712 = 0.3 Mj. 
The position of the Mode I and Mode II stationary solutions 



was recalculated for a given mass ratio and is shown by the 
thick red curves on the corresponding graphs. 

The initial conditions of the run #3 were chosen similar 
to those of the run #2: the system started in the vicinity 
of the Mode I, with the positive dissipation rate, < di > = 
10 _9 AU/yr, resulting in convergent orbits and in the loss 
of the secular energy by the system. The planetary path 
obtained is shown by solid lines in the top and middle planes 
in Figure \^-ight column. 

In analogy with the evolution of the path #2, the sys- 
tem moves in the direction of high eccentricities and of small 
ni/ri2. The smooth evolution of the system is interrupted 
by the passages through the main mean-motion resonances, 
when the planet eccentricities are strongly excited. When 
the high-eccentricity system approaches the low-order 5/1 
resonance, it is broken out. The only contrast with the path 
#2 is that, at condition y '0,1/0,2 > m,2/mi, the motion fol- 
lows closely the Mode I family, that can be clearly seen on 
the bottom graph, where time variation of the secular angle 
Azu is shown. 

The run #4 was calculated with the dissipation rate 
< di >= — 10~ 9 AU/yr and with initial conditions similar 
to those of the run #1. It is shown by black dots in Figure 
fright column. In contrast with the path #1, in this case, 
the divergent planetary eccentricities suffer strong excita- 
tions. The behaviour of the secular angle Azu on the bottom 
panel in Figure [3] shows that the system leaves quickly the 
oscillation motion around Mode I and evolves into the cir- 
culatory regime of motion. The resulting path of the run #4 
shows the capture of the system into the steady state of the 
Mode II and the consequent circularization of the planetary 
orbits. Once again, the comparison of the results of the runs 
#3 and #4 show the inadequacy of the simulation of the 
past of dissipative dynamical systems through backward in 
time integrations. 



3.3 Dependence on the dissipation rate 

In this experiment, we simulate the migrating motion of the 
planets using different values of the dissipation rate. The 
physical parameters of the system and its initial configu- 
ration are same in all simulations: the masses are fixed at 
mi = l.OMEarth and m.2 = 3.0 A/Earth, the semi-major axes 
at ai — 0.0455 AU and 02 = 0.1 AU, and the system starts in 
the vicinity of the Mode I of motion, with Azu — 0. The only 
parameter which is different in each simulation is the dissi- 
pation rate of the orbital decay of the inner planet, < di >, 
which was kept constant during each simulation. 

Figure [5] shows the planetary paths on the representa- 
tive plane (ei,e2) (top panel) and the time evolution of the 
secular angle Azu (bottom panel). The different colors cor- 
respond to the different values of the dissipation rate: cyan 
color to -3 x 10~ 8 AU/yr, red color to -3 x 10~ 9 AU/yr 
and blue color to —3 x 10~ 10 AU/yr. The thick gray curve 
represents the stationary solutions of the Mode I of motion. 
The starting configuration of all orbits is shown by a black 
symbol. 

According to our model, the migrating system will con- 
verge to and follow the Mode I family of stationary solutions, 
under the condition that the dissipation is sufficiently slow. 
This behaviour is observed in Figure [5j which shows also de- 
viations of the planet paths from the predicted route (thick 
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Figure 6. Dissipative evolution of the system along Mode I family 
of stationary solutions, for three different values of the individual 
mass of the inner planet. The mass ratio is fixed at 3, the dissipa- 
tion rate is fixed at —3 X 10 — 8 AU/yr and all orbits start at the 
same initial configuration, shown by a black symbol. The excita- 
tion due to the passages through main mean-motion resonances 
is more strong for more massive planets. The deviation from the 
conservative case is more important for less massive planets. 

gray curve). The vertical spreading of the points is origi- 
nated by the excitation of eccentricities during the passages 
through mean-motion resonances and have been discussed 
in the previous section. 

Our attention is focused on the monotonous deviations 
of the migrating paths from the Mode I family, observed at 
low eccentricities. This type of deviation is clearly related 
to the magnitude of the dissipation rate: faster is dissipa- 
tion, larger is the deviation. The final configuration of the 
rapidly migrating system consists of the circular orbit of the 
inner planet, the non-circular orbit of the outer planet (top 
panel in Figure O and the secular angle Azu, which is seen 
to be 'captured' at 90° (bottom panel). This behaviour of 
the system is associated to the phenomenon of the aperiodic 
damping, which occurs when the dissipation rate exceeds the 
proper frequency of the system (see Section iBlj) . It is worth 
observing that, if the migration is stopped after the system 
has attained its final configuration, the planets start the cir- 
culatory motion around a corresponding center belonging to 
the family of the stationary solutions. 

As will be shown in the next section, the proper fre- 
quency of the secular system decreases following a power 
law, as the eccentricities decrease (Figure [Jj. During the 
rapid decay of the inner planet, the magnitude of the proper 
frequency may become comparable to that of the dissipation 
rate. In this case, the system deviates strongly from the con- 
servative stationary solutions and can even evolve into the 
singularity at Azu — 90°, as shown in Figure [B2l in Section 
IB1I We can clearly observe these effects in the evolution of 
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Figure 7. The proper secular frequency calculated along the 
Mode I family shown in Figure [6] in the logarithmic scale. The 
mass ratio m^/mx is kept fixed at 3, but the individual planetary 
masses, initially fixed at mi = 1 AfEarth an d m 2 = 3 M^arthi were 
multiplied by factors 5, 10 and 30. The values of the mass fac- 
tor are indicated beside each curve. The vertical bands are due to 
perturbations during the passages through the main mean-motion 
resonances (MMR). 

the secular angle Aw on the bottom panel in Figure [5] It 
should be emphasizing that, in this work, we present only the 
qualitative analysis of the effects of the aperiodic damping 
on the orbits of the migrating planets. The detailed study of 
this process requires an introduction of a specific migration 
mechanism and is out of the scope of this paper. 

3.4 Dependence on the individual planetary 
masses 

In this section we perform several simulations fixing the mass 
ratio, but varying individual planetary masses. The mass 
ratio is fixed at mi/mi = 3.0, while the dissipation rate 
< di > is fixed at —3 x 10 -8 AU/yr. All simulations start at 
the same initial configuration, close to the Mode I stationary 
solution, at a\ — 0.0455 AU and 02 = 0.1 AU. According to 
our model, for given mass and semi-major axes ratios, this 
migrating system will evolve along the secular Mode I of 
motion (see Figure [2}. 

Figure [6] shows the planetary paths obtained: with the 
individual masses fixed at mi = 1 A^Earth and rri2 = 3 MEarth 
by cyan color, with these masses multiplied by factor 5 by 
red color, and multiplied by factor 30 by blue color. The 
correlation between the planet masses and the monotonous 
deviation from the Mode I family (thick gray curve) is clear: 
smaller is the mass factor, larger is the deviation of the mi- 
grating path from the route traced by the Mode I stationary 
solutions. This dependence can be easily understood from 
Figure [7] where we plot the proper frequencies obtained 
along the Mode I family. (The proper frequency is defined by 
the angular velocity of the Aro-precession.) The frequencies 
were calculated numerically for the systems starting in very 
close vicinity of the Mode I of motion (the behaviour of the 
proper frequencies around Mode II is similar). 

We note in Figure [Jj that the proper frequencies fol- 
low roughly a power law as a function of the planetary 
eccentricities. (This property allows us to understand the 
dependence on the dissipation rate described in the previ- 
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ous section.) The frequencies depend strongly on individual 
planetary masses: for small masses (~ 1 MEarth), the proper 
frequency is roughly two orders lower than for masses of sub- 
Saturn planets (~ 30 MEarth)- Thus, the large deviation of 
the Earth-like planet path from the Mode I family observed 
in Figure [5] can be explained by the effect of the aperiodic 
damping described in the previous section. 

The passages of the system through mean-motion reso- 
nances are also observed in Figure [6] In Figure they are 
associated with strong perturbation of the smooth evolution 
of the secular frequencies. It is worth noting that the per- 
turbation by mean-motion resonances is stronger in the case 
of the more massive pair of planets. 



4 THE ORBITAL ANGULAR MOMENTUM 
TRANSFER 

In the previous sections, the study was done assuming 
that the orbital angular momentum of the system is con- 
served during migration. Generally, this assumption is too 
strong and many migration mechanisms can alter the an- 
gular momentum of the system. For example, during tidal 
interactions, the rotation angular momentum of the cen- 
tral star/planet is transferred to the planet /satellite system. 
During disc-planet interactions, AM is exchanged between 
the planet and the disc, while, during the scattering of plan- 
etesimals by the planets, the angular momentum is removed 
from the system by the small bodies, etc. 

Each physical process is characterized by a specific law 
of the AM-variation. In our approach, we generalize the pro- 
cess, introducing AM-leakage as follows. Let Aai be the 
change of the semi-major axis of the migrating inner planet 
over a small time interval At. After At, the orbital angu- 
lar momentum of the system is changed by the amount 
AAM a = [AM(ai + Aai) - AM(cn)], where AM is also a 
function of the outer planet semi-major axis and eccentric- 
ities. The total variation of the orbital momentum, due to 
migration over At, is AAM = AAM a +AAM e , where AAM e 
is a contribution due to the change of the inner planet ec- 
centricity. If the angular momentum remains constant dur- 
ing migration of the system (i.e., AAM = 0), the incre- 
ment/decrement AAM a is absorbed by the system, accord- 
ing to the Equation (|10[) . and we say that there is no leakage 
of AM in the system. 

On the other hand, some portion of AAM a may 
be transferred from/to the system to/from the external 
medium (e.g., the gaseous disc, the belt of planetesimals 
or the central body rotation). If a denotes the fraction of 
the loss/gain of the orbital angular momentum, the AM- 
variation over At is 

AAM = aAAM„; (14) 
thus, the condition (|10|) must be re-written as 

Aer = (l- Q ) (1 ~ e?) Aai. (15) 

Note that, defined in this form, the AM-leakage can be pos- 
itive or negative, depending on the value of a. For instance, 
during the orbital decay of the inner planet, the system will 
lose the angular momentum, if a > 0, and gain it, if a < 0. 
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Figure 8. Evolutionary tracks given by Mode I and Mode II 
on the representative planes (ei,e2) (top panel) and (ei,ni/ri2) 
(bottom panel). Both planetary masses are equal to l.OMj and 
a2 = 1.0 AU. The leakage of the orbital angular momentum of 
the system is given by the value of a beside each curve. Arrows 
show the direction of divergent orbits. The current position of the 
system, with m2/m\ = 1.0 and orbital elements a\ = 0.22 AU, 
a2 = 1.0 AU and e\ = 0.1 and e2 = 0.2 is shown by a star symbol. 
The routes with a > are truncated at ni/ri2 < 12. 

Thus, varying a, we can use the expression (|15p to model dif- 
ferent laws of angular momentum exchange noted in planet- 
disk interactions (see Goldreich and Sari 2003, Beauge et al. 
2006) or in tidal interactions (see Rodriguez et al. 2011a). 
To illustrate the dependence on a, we assume initially that 
q is constant during the whole migration process. 

For a — 0, AM of the system is conserved during mi- 
gration; this case was studied in the previous sections. In 
Figure [8] we show the migration route composed of both 
Mode I and Mode II stationary solutions, with no leakage of 
AM (marked by 0), on the representative planes (ei,e2) and 
(ei,ni/ri2). The route was constructed for the system with 
the following parameters: mi = 7712 = 1.0 Mj, a\ = 0.22 AU, 
ai — 1.0 AU, ei = 0.1, ei = 0.2. The planets were placed on 
the aligned orbits, with Auj — 0; the position of the system 
on both graphs is shown by a red star. 

According to Equation (|14p . for a > 0, the AM is ex- 
tracted (injected) from the system during migration when 
Aai < (> 0). Four evolutionary tracks shown in Figure 
[5] were constructed with a > 0; from these, the routes with 
< a < 1.0 (0.6 and 0.9) also satisfy the condition (|13[) . In 
these cases, the divergent orbits lose the angular momentum; 
the system slides the corresponding routes in the direction 
of low eccentricities on the (ei,e2)-plane and the higher val- 
ues of ni/n2 on the (ei,ni/ri2)-plane. If we assume that, 
to attain the current position (marked by a red star), the 
divergent planets acquired one of these routes, their starting 
eccentricities must be higher when compared to the current 
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ones. It is clear from Figure [8] that, for larger values of a, 
the range of the possible starting eccentricities is smaller. 

For q = 1.0, Aei x = 0, for any Aai, and the whole 
variation of AM produced by the orbital decay is removed 
from the system. In this case, the dissipative force produces 
no change of the eccentricity, analogous to Malhotra's model 
(1995) of planet-planetesimal interactions. As described in 
Section [B2I when Ae™ = 0, the damped oscillations of ec- 
centricities and the consequent trend of the system to the 
aligned or anti-aligned configuration do not occur. Never- 
theless, as seen in Figure the guiding route of possible 
migration (marked by 1.0) still exists. 

In the case of a > 1.0, the variations Aai and Aef x 
have opposite signs, according to Equation (|15[1 . The de- 
cay (expansion) of the inner planet orbit is accompanied by 
increasing (decreasing) of eccentricities. This is because, in 
this case, the migration mechanism introduces an additional 
decrement (increment) of the angular momentum to the sys- 
tem, acting directly on ei . One such route constructed with 
q = 1.1 is shown in Figure [5] the direction of the evolution 
of the divergent system along the routes is shown by arrows. 
It should be emphasized that, in this case, the domains of 
stability of the stationary solutions are opposite to those 
presented in Figure [2] 
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Figure 9. Same as in Figure[8] except the families were calculated 
assuming all possible values of 02 and a fixed a\ = 0.02 AU. 



Finally, Equation (|15|) shows that, for a < 0, the system 
gains/loses the angular momentum even in the case of the 
orbital decay /expansion of the inner orbit. One evolutionary 
track corresponding to a = —2.0 is shown in Figure [8] The 
evolution of the system in this case is similar to the evolution 
along the routes with ^ a < 1, but the possible starting 
values of the planet eccentricities would be very high. 

Summarizing the information, we can say that, for 
a < 1.0, the migrating system, characterized by the de- 
creasing orbital energy and the increasing secular one (i.e. 
Aai < 0) , moves in the direction of the low-eccentricity do- 
main, with circularized orbits at its final configuration. On 
the contrary, when the system gains the orbital energy and 
loses the secular one (i.e. Aai > 0), it moves in the direc- 
tion of high eccentricity domains, where the close approaches 
between two planets and low-order mean-motion resonances 
destabilize the planetary motion. For a > 1.0, the direction 
of the migrating system is inverted. 

Finally, it is worth emphasizing that the migrating sce- 
nario with a constant value of a was assumed in this sec- 
tion only with the illustration propose and must be analyzed 
carefully for each specific migration mechanism, for instance, 
as has been done in (Rodriguez et al. 2011b). In that paper, 
the authors have considered three different migration mech- 
anisms, namely i) tidal interactions of the planetary system 
with the central star, ii) tidal interactions of a system of 
satellites with the central planet, and iii) gravitational in- 
teractions of the planetary system with the disc, modeled 
with a Stokes-like non-conservative force. The results ob- 
tained have shown that the AM-leakage factor a is a func- 
tion of the planet eccentricity, specific in each case. In the 
case of star-planet and disc-planet interactions, the typical 
a-values are smaller than 1. 



5 VARIATION OF THE OUTER SEMI-MAJOR 
AXIS 

Some migration processes can be described assuming that 
the dissipation affects only the outer planet orbit. For in- 
stance, this assumption is in good accordance with hydro- 
dynamical simulations, which show the outer planet driven 
inward due to torques exerted on it by outside nebular ma- 
terial (Kley 2000, 2003). Our model can be easily adopted 
to this case: The possible evolutionary tracks which the mi- 
grating system could take to attain its current configuration, 
are obtained through the calculation of Mode I and Mode 
II solutions of the Hamiltonian ((4]), for all possible values of 
a-2, and with a fixed 01. 

Figure shows by a star symbol the current position 
of a hypothetical two-planet system defined by the planet 
masses mi = U12 = 1.0Mj up , the semi-major axes ai = 0.2 
and a,2 — 0.7, the eccentricities e\ = 0.12 and e2 = 0.2, and 
the aligned periastra, Auj = 0. In analogy with what has 
been done in the case of the migrating inner planet, we cal- 
culate the evolutionary tracks parameterized by several val- 
ues of the parameter a; they are present in Figure [9] Recall 
that a characterizes the AM-leakage in the system, when 
q = corresponds to the case of the conservation of the 
orbital angular momentum of the system during migration. 
To model the AM-leakage in the system with the migrating 
outer planet, we used the Equation (|10[l . with the index 1 
replaced by the index 2. 

Figure |9j when compared to Figure |8j shows that the 
evolutionary routes of the migrating system are similar in 
both cases. Moreover, the domains of stability of the secu- 
lar Modes, I and II, which determines the mode which will 
guide the migrating system, are the same (see Figure[5}. The 
main difference is that divergent orbits are now associated 
with increasing eccentricities, if the increase of the orbital 



angular momentum due to increase of ai is not lost totally 
(a < 1.0), according to the condition (|15[) with the index 1 
replaced by 2. Therefore, the migrating systems character- 
ized by the simultaneously increasing orbital and secular en- 
ergy, move in the direction of the high eccentricity domains, 
where the close approaches between two planets could orig- 
inate a large-scale instability. The situation is inverted for 
convergent orbits, when the systems lose both the orbital 
and secular energy and, for a < 1.0, move in the direction 
of the low-eccentricity domain, with totally circularized or- 
bits at their final configurations. 

We simulated the migration of a hypothetical system, 
assuming that the outer planet experienced inward migra- 
tion losing its orbital energy. For the sake of simplicity, 
we assume that, during migration, the orbital angular mo- 
mentum of the system is conserved, i.e., a — 0. The ob- 
tained migration path is shown in Figure 1101 The system 
with mi = ni2 = 1.0 Afj up started at a\ = 0.02 AU and 
a,2 — 0.049 AU (ni/ri2 = 3.83), oscillating around Azn = 0; 
the starting position is shown by a black symbol in Figure 
1101 Evolving from the oscillation to the circulation regime, 
the system was converging gradually to the anti-aligned con- 
figuration (Mode II of motion) and finished its evolution on 
the totally circularized orbits of the planet, at n\/n,2 = 3.6. 
The time evolution of the secular angle Aw shown on the 
bottom panel in Figure [101 illustrates how the system, which 
was initially oscillating around 0, was smoothly evolving into 
the Mode II of moti on with small oscillations around 180°. 
Given that \/ aijaz < milm-\ during the evolution the sys- 
tem, this results is in good agreement with our model (see 
Figure [2]). 



6 DISCUSSIONS 

In this paper we have presented a simple method for a qual- 
itative study of the secular dynamics of the two-planet mi- 
grating system. Our model is valid to describe the evolu- 
tion of the planets which are orbiting on the same plane, 
avoiding any mean-motion resonance, and are not too close 
to produce strong short-period interactions between them. 
The evolution of the planets was studied from the point of 
view of the variation of integrals of motion of the conser- 
vative secular problem: energy, orbital angular momentum 
and planetary semi-major axes. In this formulation, no par- 
ticular dissipative mechanism was introduced; therefore, the 
results obtained are valid for any migration process. As a 
consequence, the detailed knowledge of the physical param- 
eters of the migration process and starting configurations of 
the system are unneeded in this qualitative analysis. The 
only assumption need to be done is that the dissipation pro- 
cess is weak and slow sufficiently. 

We have shown that, under this assumption, the evolu- 
tionary routes of migrating planets in phase space are traced 
by stationary solutions of the conservative secular problem. 
Parameterized by mass ratios, the families of the stationary 
solutions were constructed for continuous values of the semi- 
major axes ratio and fixed values of the angular momentum 
of the system, assuming initially that it is conserved dur- 
ing migration. We have shown that each family is composed 
of two distinct branches, corresponding to the Mode I and 
Mode II solutions, which are connected only at the origin. 
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Figure 10. Numerical simulation of the migrating two-planet 
system. Top: The starting position of the system is shown by a 
black symbol, while its final position by a red symbol. Families of 
stationary solutions are plotted by thick gray line. Arrows show 
the direction of the evolution of the system. Bottom: Time evo- 
lution of the secular angle Aw. 

During migration, the system is attracted to and follows one 
of these branches, which is said to be stable. In the case of 
divergent orbits, the Mode I (aligned orbits) is stable when 
the mass ratio and the in stantan eous semi-major axes ra- 
tio satisfy the condition y/ai/a^ < m2/mi; in the opposite 
case, the Mode II (anti-aligned orbits) is stable. For conver- 
gent orbits, the stability of the secular modes is inverted. 

We have further shown that the effects of the gain/loss 
of the orbital angular momentum of the system during mi- 
gration modifies the evolutionary tracks in the phase space. 
Knowing that the angular momentum exchange is specific 
for each kind of dissipative forces, we generalize our study, 
introducing AM-leakage. The AM-leakage is defined through 
the factor a (see Equation l |14| l). where a is a fraction of the 
angular momentum variation produced by migration (i.e. by 
the variations of Acti and/or Aa,2), which is extracted (or 
added) from the system during migration. When there is no 
AM-leakage in the migrating system (i.e. a = 0), the vari- 
ation of AM due to migration is totally absorbed by the 
system, in such a way that its orbital angular momentum 
is conserved. In this case, the eccentricity of the migrat- 
ing planet is affected according to Equation (|10|l . The case 
of a = 1 corresponds to loss/gain of the total amount of 
the AM-variation produced by migration and, in this case, 
the eccentricities of the planets are not affected by external 
forces. The values of a greater then 1 are characteristic of 
processes which produce trends of the semi-major axis and 
the eccentricity with opposite signs. It is worth noting that 
the range of theoretical values of a is large; however, the 
analysis of the orbital angular momentum exchange in sev- 
eral migrating systems driven by specified non-conservative 
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Figure 11. Schematic view of possible one-planet migration sce- 
narios in two-planet systems. 



forces (Rodriguez at al. 2011b), has shown that, in practice, 
the typical values of a belong to the interval below 1. 

Analyzing the evolutionary routes in the phase space 
we can i) conjecture on the past evolution of the system 
to its present location in the phase space; and ii) do robust 
predictions about the possible final configurations of the cur- 
rently evolving system. We have shown that the evolution 
of the secular system, under weak dissipative forces, and its 
final state depend on the balance between the orbital and 
secular energy variations of the system. To clearly illustrate 
this feature, we decompose possible migration processes into 
four elementary components; they are schematically shown 
in Figure [TT1 

The scenario 1 describes approximately the evolution 
of the system interacting with the outer protoplanetary 
disc. The planet-disc interactions are frequently modeled 
using the non-conservative Stokes force, which is charac- 
terized by the AM-leakage factor a < 1 (Rodriguez et al. 
20011b). During migration, the system loses both the or- 
bital and secular energy; as a consequence, both orbits are 
circularized during the secular evolution, while their mu- 
tual distance decreases, allowing a smooth entrance and a 
nearly zero-amplitude evolution inside a low-order mean- 
motion resonance, for instance, the 2/1 resonance. Actu- 
ally, this is the most populated resonance, with the five 
known planet pairs evolving inside it: GJ876 c-b, HD40307 
c-d, HD73526 b-c, HD82943 c-b and HD128311 b-c. An- 
alyzing these systems, we have found th at the ir parame- 
ters systematically satisfy the condition */ ai/ai < m% /mi . 
Thus, if these systems evolved according to scenario 1, they 
would converge to the Mode II of the secular motion and 
approach the 2/1 resonance in the anti-aligned configura- 
tions on low-eccentricity orbits. As shown in (Michtchenko 
et al. 2008 a,b) , this configuration is the only way to enter 
smoothly into the strong mean-motion 2/1 resonance. Sev- 
eral numerical simulations performed (Lee and Peale 2002, 
Ferraz-Mello et al. 2003, Beauge et al. 2006, Hadjidemetriou 
and Voyatzis 2010, among many others) confirm the de- 
scribed behaviour of migrating systems. 

The case 2 describes the orbital decay of a close-in 
planet due to tidal interactions with the slowly rotating 
central star. In this case, the AM-leakage factor a < 1 



(Rodriguez et al. 2011b) and the system loses the orbital 
energy and gains the secular energy; as a consequence, both 
orbits are totally circularized, while their mutual distance 
is increased. The planets b and c of the CoRoT-7 system 
have probably reached this configuration (Ferraz-Mello et 
al. 2011). The investigations of the behaviour of the hypo- 
thetical systems HD 83443 (Wu and Goldreich 2002) and 
Gliese 436 (Batygin at al. 2009) have provided the results in 
complete agreement with our model. 

An interesting example is found in our Solar System: 
the coupled motion of the longitudes of peristra of the pair 
Jupiter-Uranus, whose difference Azu — — vjj oscillates 
around 180° with an amplitude about 70° degrees (Milani 
and Nobili 1984). We can show that this peculiar configu- 
ration of the outer planets could contain a record of planet 
migration in the past. Indeed, Fernandez and Ip (1984, 1996) 
have shown that, during the scattering of planetesimals in 
the final stage of the Solar System formation, the giant plan- 
ets experienced the migration evolution: three outer giant 
planets experienced an outward displacement, while Jupiter, 
as the innermost giant planet and main ejector of bodies, mi- 
grated sunward. As shown in Milani and Nobili (1984), the 
outer Solar System can be represented by two main sub- 
systems, Sun- Jupiter-Saturn and Sun-Uranus-Neptune, ex- 
changing the angular momentum over a period of 1.1 Myr, 
same of the Aro-oscillation. In this three-body approxima- 
tion, the orbits of two subsystems were divergent during mi- 
gration and the whole system gained the secular energy. In 
the meantime, the system lost its orbital energy, due to the 
dominant mass of Jupiter. This scenario corresponds to the 
scheme 2 in Figure 1111 when the diverge nt system, whose 
parameters satisfy the condition yj a\ja% > m-z /mi , is cap- 
tured in the Mode II of motion, oscillating around 180° . 

There are not many examples of the planetary systems 
to illustrate the migration scenario 3, when the system gains 
both the orbital and secular energy. In this case, the mu- 
tual planetary distance continuously increases, together with 
the eccentricities of the planetary orbits (under condition 
a < 1). Possibly, the exosystems characterized by large mu- 
tual distances, sometimes referred to as hierarchical systems, 
have experienced this kind of evolution in the past, and this 
fact could explain consistently high eccentricities of the plan- 
etary orbits in such systems. On the other hand, the orbits 
of the hierarchical planet pairs do not exhibit the capture 
inside of one of the secular modes of motion characterized 
by oscillations of the secular angles around or 180°; thus, 
this point needs to be further addressed in more detail. 

Finally, during migration described by the scenario 4 
in Figure ITTl the system gains the orbital energy and loses 
the secular one. For a < 1, the eccentricities of two conver- 
gent orbits continuously increase during migration and the 
secular system is ultimately disrupted when it approaches 
a low-order mean-motion resonance. Thus, it seems to be 
highly improbable to find real systems which have experi- 
enced the large-scale evolution corresponding to the case 4, 
at least for the values of the AM-leakage factor a < 1. The 
picture is different in the case of a > 1, which is characteris- 
tic for tidal migrations in planet-satellite systems, specially 
those formed by giant planets (see Rodriguez at al. 2011b). 
In this case, the convergent orbits will be accompanied by 
decreasing eccentricities, similar to what occurs in the case 
1. 



Secular model of migrating planet pairs 15 



As a final consideration, we discuss the advantages and 
disadvantages of the method introduced in this paper. The 
main advantage, in our opinion, is that the method allows 
us to study the secular motion of the system evolving un- 
der a generic mechanism of migration. Thus, the patterns 
of the secular behaviour of the migrating system obtained 
in such a way, are universal. The model provides a robust 
prediction of the possible starting and final configurations of 
the systems undergone to migration processes, despite the 
qualitative approach used in its construction. The principal 
restriction of the model comes from the same fact that an 
unspecified type of dissipative forces was invoked for the dy- 
namical study. Indeed, each migration process is described 
by specific laws of the time variations of the energy, the 
orbital angular momentum and the semi-major axes of the 
planets. Treating the problem in general terms of variation 
of these quantities, we lose the information on the migration 
time-scales, characteristic of the system under study. Thus, 
to obtain a complete solution of the problem, our model 
must be completed by additional analytical (e.g. Mardling 
2007) or numerical (e.g. Rodriguez et al. 2011) investiga- 
tions of the time-dependence of the secular evolution on a 
specific dissipative force. 

APPENDIX A: THE BASIC CONCEPTS OF 
THE CONSERVATIVE SECULAR DYNAMICS 

Al Secular model 

In order to understand the secular behaviour of the planet 
pair, we study the second-order (i.e., 0(e 2 )) approximation 
of the Hamiltonian describing the secular planar three-body 
system reduced to one degree of freedom (for details see 
Michtchenko & Ferraz-Mello 2001, Callegari et al. 2004). 
Discarding constant Keplerian terms, the reduced secular 
Hamiltonian can be written, in terms of the canonical elliptic 
variables (J3J), as 

"Hscc = 2(C - D)h + 2£^7i(AMD - h) cos Aw, (Al) 

where AMD is given in ((5]) and the coefficients are defined 
as 

C = -(2)(°) Ms/L^l 

D = -(3)(°»M 5 /I|, (A2) 

e = -(2i) ( - i) M 3 /Li^rrx;. 

The coefficients (2) (0) = (3) (0) and (21) (_1) are functions 
of the Laplace coefficients bi k \ai/a2) written in the form 
adopted by Le Verrier (1855) (for their explicit expressions, 
see Callegari et al. 2004) and Mz = fJ% m i m 2 ■ 

The long-term dynamics of two massive planets evolv- 
ing far from any mean-motion resonance is defined by the 
averaged Hamiltonian ()A1|) ; in the non-dissipative approach, 
■Hsec is the secular energy of the system, which is conserved 
along the motion. As the energy Usee, the orbital angular 
momentum ((6J and the angular momentum deficit ([5]) are 
also constants of the of the secular planet motion (they were 
used to reduce the usual two-degrees-of-freedom Hamilto- 
nian) . 

The AMD appearing as a parameter in the expression 
of the secular Hamiltonian (|A1|) has a clear algebraical inter- 
pretation: it has a minimum value (zero) for circular orbits 



and increases with increasing eccentricities. The behaviour 
of AM is reverse of that of AMD. It should be emphasized 
that the AMD is conserved only in the secular problem be- 
cause, in this problem, L\ and 7/2 (i.e., a± and 02) are con- 
stants of motion. Also, the action variable 7i and its anal- 
ogous I2 given in ((3} define the partial angular momentum 
deficits of the planets (their sum is the AMD). 



A2 Phase space 

In order to have a geometrical representation of the phase 
space of the secular system, we transform the Hamiltonian 
(| All) following the approach introduced by Pauwels (1983). 
First, we note that the coefficient E given in (|A2I) is always 
positive, while C and D axe always negative. We introduce a 
new variable 8 and new parameters a and /?, in the following 
way: 

P cos a sin a 

"Hscc = 2 (C^~D) h +2 "iT y/h (AMD - h) cos Avu, (A3) 
v -v' v ' 

AMD(l+COsS) AMD sin 6 

where both a and 8 vary in the range from to n. We also 
have P 2 = {C-D) 2 +E 2 and Equation (|A"3| then is rewritten 
as 

Hscc-P AMD cos a = ft AMD(cosqcos <5+sin a sin <5 cos Aw).(A4) 

The variables (<$, Aw) may be considered as defining a 
system of spherical coordinates, with a vertical polar axis 
NS, on a sphere (Pauwels' sphere) of radius ft AMD. The 
right-hand side of Equation (|A4[) is the expression of a cosine 
law and can be re-written as 

£1 = Usee - ft AMD cos a = ft AMD cos 8* , (A5) 

where S* is the distance of the generic point P with spherical 
coordinates (8, Avu) to the pole Z, defined by the axis Z Z' , 
which is inclined of a with respect to the vertical N S (see 
Figure [AT]) . The points with 5* = const form parallel circles 
on the sphere, whose centers belong to the axis Z Z' . The 
intersections of the axis Z Z' with the sphere have spherical 
coordinates 8 — a and Au? = 0, for Z-pole; for Z'-pole, 
8 — 7r — a and Azu = tt. These intersections correspond 
to the two equilibria of the one-degree-of-freedom secular 
Hamiltonian ()A5|) . 

8* is the polar angle in a system of spherical coordinates 
whose axis is Z Z' (instead of the vertical axis). In analogy 
with the transformation given in Equation (|A3|) . we may 
introduce a new variable /* through 

I* = ^^(l + coscT) (A6) 

and write the secular Hamiltonian (|A5|I as 

£1= p (27* — AMD). (A7) 

Some important consequences of the performed trans- 
formations come from the fact that Jj* may be considered as 
a new action variable. The conjugate of 7* is an angle wi, 
rotating with the constant angular velocity u>i = 2/3, which 
is independent on the energy. As a consequence, the periods 
of the motions on all parallel circles are same and depend 
only on /3, i.e., on the constants C, D, E. 
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Figure Al. Separation of the Ara modes of motion on Pauwcls' 
sphere. In the spherical caps defined by parallel circles passing 
through the N and S poles of the sphere, Aro oscillates: around 
in the cap with Z vertex and around n in the cap with Z' vertex. 
In the zone between the two parallels, Aro circulates. Left: Case 
< a < it/2. Right: Case n/2 < a < it. 




Figure A2. Locus of the curves I* = const in the plane 
II cos Aro, 7i sin Aoj. The parallel circles shown in Figure fXT] ap- 
pear here as continuous lines. In both cases, I* is maximum at the 
Mode I center. Left: Case < a < n/2. Right: Case n/2 < a < n. 



The new action variable II (as well as the polar an- 
gle 6*) is a constant of the motion. Since (3 AMD > 0, the 
secular energy E\ and, therefore, H*scc is maximal (resp. 
minimal ) at the pole Z, where 5* = (resp. Z' , where 
5" = n). 

The dynamical interpretation of the solutions defined by 
the secular Hamiltonian is now simple. For a given AMD, 
the solutions are the curves on the sphere corresponding to 
points of the same energy, that is, parallel circles. Indeed, 
the conservation of the energy and of the orbital angular 
momentum (or the angular momentum deficit) along one 
solution define the main features of the secular motion of 
planetary systems. 

The solutions in the spherical zone between the two 
parallel circles, defined by 5* — a and 8* — n — a (see 
Figure |A"T|) . circulate around the vertical axis N S and, along 
them, the secular angle Azu circulate from to 2n. In the 
solutions inside the spherical caps between these parallels 
and the poles Z and Z' , the angles Azu remain bounded. 
In the spherical cap whose vertex is at Z (resp. Z'), the 
angle Azu remains in the interval —n/2 < Azu < n/2 (resp. 
n/2 < Azu < 3n/2) the perihelia move in such a way that 
Azu oscillates around (resp. n). It is worth noting that 
this behaviour does not depend on which hemisphere on the 
Pauwels' sphere each pole is located. 



For the negative coefficients C and D, when \C\ < \D\ 
(resp. jCj > \D\), a < n/2 (resp. a > n/2) and the pole 
Z appears in the upper (resp. lower) hemisphere of the 
Pauwels' sphere (see Figure IaT|) . The transition between the 
two hemispheres occurs when C = D (or L\ — 7/2 ) and, in 
this particular case, both centers are located on the horizon- 
tal plane of the Pauwels' sphere. This last condition can be 
rewritten, up to second order in masses, as in (|12|) . 

Finally, let it be said that the derivations given above 
can be done using the variables (—I2, Azu). Indeed, it is 
sometimes useful to use simultaneously the two possibili- 
ties for better characterizing the planets dynamics. The re- 
sults are equivalent to the above ones and are obtained just 
putting I\ = AMD — I2 in the equations. 



A3 Stationary solutions 

As shown in Figures [A 1 1 and TA2I two possible mutual orien- 
tations of the planet orbits are of special interest to under- 
stand the secular dynamics: one is the case of aligned orbits, 
when Azu = 0, and the other is the case of anti-aligned or- 
bits, when Azu = 180. These particular configurations cor- 
respond to the stationary solutions of the secular three-body 
(two planet) problem and have been dubbed as Mode I and 
Mode II (Tittemore and Wisdom 1988, Michtchenko and 
Ferraz-Mello 2001). 

The Mode I corresponds to the center located at the 
pole Z on the Pauwels' sphere, while the Mode II corre- 
sponds to the center located at the pole Z' . It is worth 
mentioning that, by the definition (|A6|I . 7* is equal to zero 
(minimum of /*) at the center II and equal to AMD (maxi- 
mum of /*) at the center I. In the alternative representation, 
when I2 is used instead of 7i , the picture is similar to Figure 
IA2I just by commuting the left and right panels. We intro- 
duce I2 = AMD — /* and then we have that 7| is equal 
to zero (minimum of 7|) at the center I and equal to AMD 
(maximum of 7J) at the center II. 



A4 The general secular problem 

It is worth stressing the fact that the secular problem defined 
by the Hamiltonian (|A1[) is a low-order approximation valid 
only for small eccentricities and large separations between 
planets. To expand the model to high eccentricities, we use 
a semi-analytical approach employing a numerical averag- 
ing of the short-period gravitational interactions of the two 
planets (Michtchenko and Malhotra 2004). The Hamiltonian 
describing precisely the secular perturbations between the 
planets on high-eccentricity orbits, is given by 

'Hscc = - 7 ^ ? / / 1 m2 R(ch, e u Mi, m) rfMidM 2 , (A8) 

O) Jo Jo ffl 2 

where the integrand is the direct part of the disturbing func- 
tion written in terms of the canonical astrocentric semi- 
major axes and eccentricities of the planets (a,i and a, re- 
spectively) and of the angular elements: mean anomalies Mi 
and longitudes of pericenter zui. 

The phase space is no longer the same as shown in Fig- 
ures I All and IA21 since those structures are now perturbed 
by the high-order eccentricity terms. The role played by the 
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Figure A3. Families of stationary solutions of the conservative 
secular problem, parameterized by the planet mass ratio (top 
panel) and the semi-major axes ratio (bottom panel). Mode I 
solutions (Aro = 0) are in the right half-plane, with positive val- 
ues on the el-axis. Mode II solutions (Aro = 180°) are in the left 
half-plane, with negative values on the el-axis. The branches of 
the curves in black color corresponds to solutions with /f < 7f 
and in red color to solutions with /f > if ■ Along each curve, 
the values of the parameters m2/m\ and a\/a2 are fixed, but the 
values of the energy and the angular momentum vary. 

cases |C| < \D\ and |C| > \D\ is now more involving. How- 
ever, the important characteristics continues to be the fact 
that the pole Z (associated with the maximum of the sec- 
ular energy) is located either in the upper or in the lower 
hemisphere. The transition between two cases occurs when 
the pole Z appears on the horizontal plane of the Pauwels' 
sphere, i.e., when 8z = tt/2 or if — |AMD. We may still 
write this condition as if = if or, up to second order in 
masses, as: 

mi _ 1- yT-ef 2 

m 2 Y a 2 i _ ^i-ef 2' 

where ef and ef are the eccentricities of the stationary solu- 
tions. The linear approximation (|12[) results from the above 
expression when ef = ef . It should be stressed, that the 
condition ()A9|) defines the mode of motion to which the sec- 
ular system evolves when dissipative effects are included. 

Figure DOl shows the families of stationary solutions pa- 
rameterized by the ratio of the planet masses and by the 
semi-major axes ratio. In the top panel, the ratio of the semi- 
major axes is fixed at ai/di = 0.2; in the bottom panel, the 
mass ratio is fixed at mi/ mi = 1. Along each curve, the 
values of both parameters are fixed, but the values of the 
energy and the angular momentum deficit vary. We calcu- 
lated the values of the partial angular momentum deficit of 
the planets for each stationary solution shown in Figure lAlfl 
The relation between them varies along each curve: when 
if < 1% 1 we plot the branches of the families in black color, 
and, when if > if , we plot them in red color. 



APPENDIX B: ELEMENTS OF THE 
DISSIPATIVE SECULAR DYNAMICS 

In this section we illustrate the main features of the dynam- 
ics of secular systems with dissipation. For this task we use 
the simple model developed in the previous section. In the 
conservative dynamics, the secular energy given by Hamil- 
tonian (|A7|I is a function of the action variable /j* (conse- 
quently, of ei). The dissipation of the energy occurs when 
friction forces are introduced and the system loses (or gains) 
the secular energy during the evolution. 

Moreover, the Hamiltonian ()A7[) is parameterized by 
the angular momentum AM (through AMD = Li+L 2 — AM) 
and the ratios m2/mi and axja^. All these parameters may 
vary under action of some external processes (e. g. planetary 
migration, mass-loss in close-in planets etc). In particular, 
the variations of the planetary semi-major axes are generally 
related to changes of the orbital energy and orbital migra- 
tion of the planets. However, it is worth noting that, in this 
section, we do not investigate the dissipation of the orbital 
energy, defined by the Keplerian motion of the planet. We 
concentrate our attention on the secular component of the 
energy given by the Hamiltonian Q , which contains no Ke- 
plerian terms (these terms being constants in the secular 
conservative problem were omitted). Nevertheless the secu- 
lar Hamiltonian is still dependent on the semi-major axes 
of the planets (more precisely, on their ratio) as a parame- 
ter; thus its change will provoke the changes in the secular 
dynamics. 

To facilitate our understanding, we vary the described 
above constants separately and study the behaviour of the 
system in each case. 

Bl Secular motion with friction 

First, we consider the case of the loss/gain of the secular 
energy in the system defined by the Hamiltonian ()A7[) and 
moving under action of a friction force. The equations of 
motion may be written as: 

H = = (Bl) 

aw* 

w = !|=2/3, (B2) 

where we assume a general form of the w*-independent 
non-Hamiltonian term /(/f ), with a coefficient of friction e. 
The function /(/*) is so-called Rayleigh's dissipation func- 
tion, given by a positive definite quadratic form in the time 
derivative of the coordinate, and represent frictional forces 
which are proportional to velocities (for details see Mierovich 
1970, Chapt. 2). The coefficient e has either positive values, 
in the case of the loss of energy in the system, or negative 
values, when the system gains energy. It is worth emphasiz- 
ing that the friction force acts on the variable 7f (ei) alone, 
while the parameters (the angular momentum, its deficit 
and the mass and semi-major axes ratios) remain constant 
during the evolution. Since AMD and 2/3 are constant, the 
motion of the damped system occurs on the Pauwels' sphere 
of the unchanged radius. 

As shown in Section lA"2l the secular system is linear on a 
sphere; thus its behaviour is similar to well known behaviour 
of the harmonic oscillator moving under acting friction force 
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Figure Bl. Top: Schematic presentation of a damped oscillatory 
process, which takes place when the dissipation rate is slower than 
the proper oscillation period. The spiral-like trajectory of the sec- 
ular system is winding to Mode I of motion (left graph), which 
plays role of the stable focus. On the right graph, the other side of 
the Pauwels' sphere shows the same trajectory, unwinding from 
Mode II, which is an unstable focus. The dashed curves are cir- 
cular trajectories of the conservative system shown in Figure fATl 
Bottom: the same trajectory shown in the top graphs, on the 
(ei ,e2)— plane. The gray thick lines are the families of conserva- 
tive stationary solutions, Mode I and Mode II. The system starts 
at the configuration show by a black symbol and is stopped at 
configuration shown by a red symbol. Arrows show the direction 
of the drift of the dissipative motion. 



(see for details Andronov et al. 1966, Chapt. I). For e = 0, 
the energy is conserved and the solution of the system (|B1|) - 
(|B2|) is trivial: I* = const and the angle w* is obtained 
through a simple quadrature as w* = 2/3 1 + const. 

If friction is positive (e > 0), I* is decreasing to its mini- 
mal (zero) value. If dissipation rate is slower than the proper 
period of the system, the decrease of I* provokes damped 
oscillations of the eccentricity e\. All possible trajectories on 
the Pauwels' sphere, which were simple circles in the conser- 
vative case (see Figure |A"T|) . have now the form of a spiral; 
one of these trajectories calculated for the system (|B1|) is 
shown in Figure IBll top. All together, they form a family of 
spirals, enclosed in each other, with an asymptotic point at 
the equilibrium of the system (|B1[I . A singular point of this 
kind is known as a focus. 

The stability of a focus is related to whether the tra- 
jectories in its close vicinity are winding (incoming) or un- 
winding (outgoing) with respect to the direction of motion. 
If the orbits are incoming, the focus is stable. For example, 
on the Pauwels' sphere in Figure [BTI /eft. the stable focus is 



Figure B2. Stationary solutions of the dissipative system (I B II) - 
||B2| | on the plane (Aro,<5). Large dots show location of the so- 
lutions of the conservative problem when e = 0. The left half- 
plane is characterized by the negative e-values, while the right 
half-plane by the positive e-values. Arrows show direction of in- 
creasing e (in abs. value). Black curve was obtained for the mass 
ratio m2/m± = 2.0, while the red curve for m2/m\ = 0.3. 



close to the Mode I of motion. Otherwise, moving along a 
spiral orbit, the system moves away from the Mode II (Fig- 
ure \Bl\ naht). in direction of Mode I; in this case, Mode II 
is an unstable focus. The same trajectory projected on the 
(ei ,C2)-plane is shown in the bottom panel in Figure lBll The 
evolution of the system on this plane is confined to the level 
of the constant AMD: it starts in the vicinity of the Mode 
II stationary solution (black cross) and is stopped when the 
system reaches the Mode I stable configuration (red cross). 

In the case of the negative friction, the system gains 
energy. There will be no longer the damping but the rein- 
forcement of the oscillations. The variable /* is increasing 
to its maximal (AMD) value. Portrait on the phase space is 
still a family of spirals, however, the direction of the spirals 
is inverted with respect to that obtained for systems with a 
positive friction. 

It is worth noting that, in contrast with the dynamics of 
the harmonic oscillator, the increasing with time deviation 
of the secular system from the unstable focus (Mode II in 
Figure IBTJ) has an upper limit, given by the position of Mode 
I of motion. This is due to the fact that the phase space of 
secular systems is a sphere; thus all possible trajectories, 
unwinding the unstable focus, will converge to the stable 
one. As a consequence, the dissipative system will ultimately 
evolve into a steady state. 

Finally, if friction is positive (e > 0) and its rate is faster 
then the proper period of the system, we have an aperiodic 
damping. This case appears in our simulations and will be 
discussed in Section 5. 

To illustrate the dependence on the friction coefficient, 
we calculate the equilibria of the system (|Bl[l - (|B2p as a 
function of e. For this task, we need a specific form of the 
friction force: for the sake of simplicity, it was chosen to be 
proportional to the action variable J*. The coefficient e was 
varied in the range from to 1, in the units of the proper 
frequency of motion, 2 /3. 

The solutions are shown in Figure lB2l in the plane (S, 
A-cu), where Aw is the characteristic angle of the secular 
model and 5 is a function of the partial angular momentum 
deficit of the inner planet, I\, and, consequently, of the inner 
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Figure B3. Trajectory of the secular system with the inner 
planet in the linear regime, on the (ei,ea)— plane. The gray lines 
show the families of Mode I and Mode II stationary solutions. 
The system starts at the configuration show by a black symbol 
and is stopped when both orbits are circularized. Arrow shows 
the direction of the drift of the center. The orbit was calculated 
with the mass ratio m,2/m\ = 2.0 and the fixed value of AM. 



eccentricity. The parameter 01/02 of the secular model was 
fixed at 0.2 and two values of the parameter milm\ were 
used: 2.0 and 0.3. For a given 01/02 and small eccentricities, 
the first value of mij 'mi satisfies the condition \C\ > \D\ 
and the second the condition |C| < \D\ (see Equation (I12|) ). 
The critical value of the mass ratio, obtained for \C\ = \D\ 
and ai/a2 = 0.2, is mijmx — 0.447. 

The black curve in Figure IB2I shows the foci obtained 
for m2/m\ — 2.0. Since |C| > \D\ in this case, the conser- 
vative equilibrium Mode I is in the lower hemisphere of the 
Pauwels' sphere, where 5 > n/2, and the opposite Mode II is 
in the upper hemisphere, where 5 < n/2. Their positions are 
shown by large dots on the black curve, at Aro = 0/360°, 
for Mode I, and, at Auj = 180°, for Mode II. 

The conservative stationary solutions (at e = 0) are ex- 
trema of the functions 5(e) and Avj(e). Therefore, for small 
non-zero values of e, the equilibria of the dissipative system 
are close to the equilibria of the conservative one. For this 
reason, for small e, we continue to refer the foci close to 
Azu = as Mode I and close to Auj = 180° as Mode II of 
motion. 

Discontinuities which occur at Azu = 90°/270° and 
S = 0/180° clearly separate two modes of motion. From 
the definitions in Equation (|A3|I . it is easy to verify that 
theses discontinuities are related to the singularities of the 
secular problem (1A1[) . which take place at Ii — I2 — and 
Ji = I2 (see Equations (10) and (11) in Michtchenko and 
Ferraz-Mello 2001). 

The foci obtained for mj/mi = 0.3 are shown in Figure 
IB2l bv the red curve. In this case, |C| < \D\ and the location 
of the equilibria is opposite to that of the previous case: the 
Mode I is located now in the upper hemisphere (8 < n/2), 
while Mode II is in the lower hemisphere (S > ir/2). In 
the rest, the evolution of the foci with the increasing value 
of e is similar; particulary, the discontinuities also occur at 
Ii = I2 = and h — 12- 

Which of the Mode I and Mode II foci will be stable or 
unstable depends on the sign of e and the parameters of the 
system and is discussed in Section TB3I 



B2 Variation of the semi-major axes ratio and the 
angular momentum 

Now we consider the dissipative system shown in Figure [Bll 
whose semi-major axes ratio varies, due to some hypotheti- 
cal processes which keep unchanged the other two parame- 
ters, 1712/ mi and AM. Although AMD is not longer constant, 
the phase space of the system can be still considered as the 
Pauwels' sphere, whose radius varies during the evolution, 
following the changes of AMD. In such a way, we introduce 
an 'instantaneous' Pouwels' sphere. 

If the ai /(^-variation is sufficiently slow, the trajec- 
tory of the system on the sphere continues to be a spiral 
approaching a stable center, but the position of the center 
on the Pauwels' sphere is continuously dislocated together 
with changes of 01/02- Projected on the (ei,e2)-plane, the 
trajectory of the system is not longer confined to one AMD- 
level (as in Figure \BT\ bottom). As shown in Figure [B3l the 
damped oscillations of the eccentricities occur now around a 
center that drifts along the family of Mode I stationary solu- 
tions of the conservative case. This family shown by a gray 
line was calculated using the Hamiltonian ([4]) of the general 
model, assuming continuous values of the semi-major axes 
ratio and the fixed values of the angular momentum and the 
mass ratio (m2/mi = 2.0). 

The trajectory shown in Figure IB3I was obtained 
through numerical integration of the equations of dissipa- 
tive motion of the secular system (|B1[) - (|B2I) . including the 
linearly decreasing in time ai . When ai is decreased (AMD 
is also decreased assuming AM =const), the center of the 
oscillating system slides down the curve of stationary solu- 
tions, in direction of the origin. The evolution is stopped 
when the system approaches the origin, when both plan- 
etary orbits are circularized. At the origin, AMD = and 
the phase space of the secular system is degenerated in a sin- 
gle point (radius of the Pauwels' sphere tends to zero). As 
a consequence, we can calculate the minimal possible value 
of 01/02 during the orbital decay (Rodriguez et al. 2011a). 
Indeed, imposing ei = e2 = in Equation (|6]), we obtain, 
up to first order in masses, 

M°f) =(AM'-l)^, (B3) 
V Va2/min mi 

where AM' = AM/n^^/aJ. 

When the orbit of the inner planet is expanded, the 
center slides up the Mode I family, in direction of high ec- 
centricities. In the end, the trajectory will reach a very-high- 
eccentricity region, where resonant and short-period inter- 
actions destabilize the planetary motion. 

The behaviour of the dissipative secular system in the 
case when the angular momentum AM varies, but the other 
parameters are conserved, is essentially same as described in 
the previous case (see Figure [B3| . The changes of AM pro- 
vokes dissipation of the energy and the angular momentum 
deficit. As a result, the centers of the secular motion drift 
on the Pauwels' sphere during the evolution of the system, 
meanwhile the sphere changes its size. The only difference 
with the previous case is that the evolutionary track of the 
system is calculated for all possible values of the angular 
momentum and the fixed 01/02 and m2/mi. 



20 T.A. Michtchenko and A. Rodriguez 



B3 Stability of the centers 

The determination of stability of the centers is immediate 
when the problem is written in terms of the variable It 
defined in Equation <|A6[) . Indeed, in this case, I{ (and the 
secular energy) is maximal (minimal) at Mode I (resp. Mode 
II) of motion (see Section IA2[) . Thus, for positive friction 
when the energy is lost, the system will converge to the 
Mode II, and, for negative friction, when the system gains 
energy, to the Mode I. 

However, the determination of stability of the centers 
becomes more complicated when the dissipation is given in 
terms of eccentricity ei, as in Equation (110p . The rigorous 
way to obtain the solution in this case is to assert the re- 
lationship between the variations IS.Il an< A Aei, performing 
inverse transformation from the variable ij* in Equation (|A6[) 
to the eccentricity e\. 

The other way to solve the problem is based on the 
definitions done in Section [X] For instance, from Equations 
(|10[) and (1A3[) . we deduce that the orbital decay Aa± < 
will produce the following variations: 

Aef < => Ah < => AS > 0. 

The above conditions show that, during the orbital decay, 
the system will evolve into the center with the smaller partial 
angular momentum deficit if , which is located in the lower 
hemisphere of the instantaneous Pauwels' sphere. On con- 
trary, during the expansion of the orbit of the inner planet 
(Aai > 0), the system will evolve into the center with the 
larger if . It is easily verified that the results are equivalent 
substituting index 1 by the index 2 in the above deduction. 
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